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Zonal  Harmonics  in  Low  Frequency  Terrestrial 
Radio  Wave  Propagation 


J.    Ralph  Johler 


If  7?n    and  Tn    are  ground  and  ionosphere  reflection  coefficient 
matrices  and   1    is  the  unit  matrix,    then  for  the  terrestrial  waveguide 

boundaries,    a  geometric  series     \l   +    2     7?^  T^  |    for  each  zonal  harmonic 

i  =o 
mode,    n,    results  in  a  rapidly  converging  zonal  harmonic   series.      Thus, 

by  an  interchange  of  the  order  of  summation  of  the  two  series,    the  zonal 

harmonic  series  can  be  employed  to  calculate  each  term  of  the  geometric 

serie  s. 

This  method  of  analysis  provides  a  practical  solution  to  the  propa- 
gated LF,    VLF,    ELF  terrestrial  waveguide  field.      Thus,    geometric 
series  terms  (often  called  wavehops)  can  be  calculated  quickly  and 
accurately  with  the  zonal  harmonic  series.      The  summation  of  the 
geometric   series  then  provides  an  efficient  method  for  calculating  the 
total  field  (or  waveguide  mode   sum).      The  simplicity  of  the  analysis  is 
an  attractive  feature  when  compared  with  the  complex  integral  method. 
Computation  speed  is  obtained  on  the  large   scale  computer  by  use  of 
recursion  formulas  for  the  Legendre  and  spherical  wave  functions.     It 
is  noted  that  the  terms  of  the  geometric  series  in  the  high  frequency 
limit  give  the  geometric -optical  rays. 

An  improved  definition  of  the  reflection  coefficient  matrix  is 
presented  which  accounts  for  the  reflection  process  in  a  manner  which 
can  be  justified  both  mathematically  and  physically.      In  fact,    the  entire 
reflection  coefficient  matrix  for  the  anisotropic  ionosphere  can  be 
retained  as  a  variable  of  the  zonal  harmonic  summation  process  with 
resultant  improvement  in  computation  accuracy.      Demonstrative  compu- 
tations indicate  considerable  advantage  in  the  method  as  an  alternate 
approach  to  the  propagation  problem. 


Key  Words:      Extra  low  frequencies,    geometric -optics,    geometric- 
series,    LF,    VLF,    ELF  mode  theory,    low  frequencies, 
terrestrial  radio  wave  propagation,    very  low  frequencies, 
zonal  harmonics. 


1 .      Introduction 

In  a  previous  paper  [johler,    1964]  the  notion  of  a  zonal  harmonic 
series  for  each  term  of  the  geometric   series  representation  of  the 
terrestrial  radio  wave  field  was  introduced.      The  upper  boundary  of  the 
terrestrial  radio  waveguide  for  low  frequency  radio  waves    (Johler 
and  Berry,     1964)     can  be  characterized  by  the  reflection  coefficient 
matrix   T, 


T  = 


T(E)  T(s) 

-1-  e  e  ■*■  e  m 

T(S)  T(s) 

in  6  in  ni 


(i.i) 


1  (s  ) 

where       Te  e    is  the   reflection  coefficient  for  vertical  polarization 
(TM-mode)  of  the  incident  and  reflected  waves  on  the  boundary  g  -  a  +  h, 
figure  1,    where  a  is  the   radius  of  the  terrestrial  sphere,    and  h  is  the 

(s) 

height  of  the  guide.      Similarly,    Tm      refers  to  the  vertical  magnetic 

(S)  (s) 

incident  and  reflected  wave s  (TE-mode).      Also,    Tem   and  Tm  e   are  abnormal 

components  resulting  from  the  anisotropic  nature  of  the  ionosphere.      For 

simplicity,    it  is  instructive  to  assume  the  ionosphere  to  be  isotropic, 

(£) 
and  the  excitation  TM.      Then,    T  =  T^  e\      Also,    for  the  ground, 


7?  = 


(s) 

r:  '         o 


R 


(■) 


(1.2) 


(s)  /s\ 

where  Re      refers  to  vertical  polarization  of  electric  vector  and  R^ 
refers  to  vertical  polarization  of  the  magnetic  vector.      For  the  TM  exci 
tation  of  the  terrestrial  waveguide  with  isotropic  reflection  at  the  iono- 
sphere, 


7?  =  R 


(s) 


(1.3) 


or  dropping  the   subscript,    e, 


i?  =  R 


(■) 


(1.4) 


Superscript  (s)  is  here  employed  to  remind  the  reader  that  the   spherT 
cal  convergence -divergence  or  focusing  factor  has  not  as  yet  been 
removed  from  the  reflection  coefficient.      See  (1.  10). 


The  classical  zonal  harmonic  solution  [johler,    1964],    can  be 
written  for  terminals  S  and  O  in  figure   1,    on  the   ground,    r  =  a,    for  the 
particular  case  of  the  vertical  electric  polarization,    Er  , 


E.   =  -^ 


'r     ~    i^2  _  4 


^     £    n(n+l)(2n+l)Pn(cos  0)  cfa    K 


kfa4     4n 

n  =0 

x{l   +   R(S>^}{1  +  Tr>^}{l-Ri!»T<S)r-  I1'5' 

bl  a 

(n  =  0,    1,    2,    3.  .  .  ) 

Here  I0  t  is  the  intensity  of  the  electrical  point  source  dipole,    S,    figure  1, 
in  ampere -meters  which  can  be  for  convenience  taken  as  unit,    10^=  1. 
The  permeability  of  space  is  |J.0,    |J.0    =    4tt(10~    )  henry /meter  and  c  is 
the  speed  of  light  (c  ~  2.  997925  (108)  m/s).      kx   is  the  wave  number, 

UL)  . 

kx    =  —  n-L    in  air  where  r^   ~   1.  0001   to   1.  0003  and  the  frequency  f  =  uj/2tt. 

uu 
In  vacuum,    ^    =1,    k0    =  —  .      The   subscript  n  on  Tn    and  Rn    reminds  the 

reader  that  these  quantities  are  calculated  for  each  n.      Pn(z)  is  the 
solution  of  Legendre's  differential  equation, 


<L    -''l^-^^^HUP.I^O.  (1.6) 

(l  ,  2  ) 

The  abbreviations  Q1&'        and  ^  a   are,    letting  kx  a  -  z, 


ka  =    Cn  (z)=^/—    Hn  +  i      (z)  I1-7) 


'ia     =    K(Z)    =      /^    Jn  +  i(z)  (1-8) 


(l      2  ) 

where  Jn+i(z)  and  H^  +  'i  (z)  are  Bessel  and  Hankel  functions  of  order 

9  2 

(12) 

nf-g-  and  argument  z  and  H^  may  be  of  the   first  or   second  kind. 
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The  reflection  coefficients  R;  '  and  T^      contain  focusing  or 
convergence -divergence  factors.      Thus,    the  spherical  reflection  coeffic 
ients  will  be  designated  by  Rn    and  Tn  ,    dropping  the   superscript  (s),    or, 
after   substitution  in  (1.5)  of  the  exact  relationship, 


*n(«)  =i[cil)(z)  +  daW] 


(1.9) 


R(S)    _     -Cia      R 
Rn        -  (a)      R> 

bl  a 


(1.10) 


T 


(2) 

(■>  _  ^cii  T 

n  "  r(l)    n 


(1.11) 


Here,    the  spherical  reflection  coefficients  are  defined, 


and, 


Rn     = 


T. 


c(1)/  k    c(1)/ 

'si  a         Ki      ^3» 


Ci 


(1) 


k2     r(i) 


c(3)'  k  c(l)'  ' 

'si  a  %_  'oSa 

r(2)  k2        (i) 

Wia  "=>2a 


ill.^ 


^3e 


(2)         k3     r(a) 


£i  I ksLa 


-c(1)/    k    c(2) 

r(i)         k3    r(2) 


36 


(1.12) 


(1.13) 


where   C~  '     '  (z)   =  f-    C1*  '  ~ '(z). 


dz 
The  wave  number,    k3 ,    is, 


k3   = 


(!) 


UU. 


uo(v+iO)) 


(1.14) 
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2  Ne2 

for  a  plasma  frequency  squared  IUN where   N  is  the  number  density 

e0m 


of  the  plasma  (N  =  Neel/m3   for  an  electron  plasma  and  m  =  me    the  elec 

tronic  mass)  and  £0    =  — .      The  electron-neutral  collision  frequency 

c    Ho 

is  v,    sec-1.      The  wave  number  of  the   ground,    k2 ,    is  given  by, 


W     /  •  °"M-  o  c2 


where  the  permittivity,    €  =   €0  £2  .      The   conductivity  is  given  by  O, 
mhos  /m. 

The  geometric  series  expansion  was  obtained  from  the  denominator 
in  (1.  5), 


i-kW 


_1=  1  +  £    [r1S)t1S)]   ,  (1.16) 


3=1 


l      (s)      (s)i 
which  converges  absolutely,    if    |  R*      T^         <  1. 

This  basic  idea  has  been  employed  in  integral  form  in  the   complex 
n-plane  by  Bremmer  [1949]  and  Wait  [  1961  ].     Indeed,    Wait  [l96l] 
suggested  a  rigorous  type  of  geometric -optics  applicable  to  low  fre- 
quency propagation.      Berry  [  1964a,  b]  and  Berry  and  Chrisman  [1965] 
developed  elegant  computation  methods  based  on  this  idea  using 
asymptotic  formulas  for  the   spherical  wave  functions  (1.  6),    (1.  7),    and 
(1.  8)  in  the  complex  n-plane.      In  effect,    the  method  presented  in  this 
paper  extends  the  classical  geometric-optics    [Johler,    1961;    1962] 
to  great  distance  (>  3000  km).      Such  an  approach  is  an  alternate  to  the 
waveguide  mode   solution  described  in  detail  by  Wait  [1962],      The 
asymptotic  formulas  for  the    spherical  wave  functions    [Berry,    1964a,  bj 
are  intricate  to  program  for  a  large   scale  computer.      On  the  other  hand, 
along  the  real  axis  of  the  complex  n-plane,    the   spherical  wave  functions 
(1.  6),    (1.  7),    and  (1.  8),    for  integer  values  of  n,    can  be   calculated  with 
great  speed  and  simplicity  on  a  computer  with  standard  recursion  formu- 
las.     Furthermore,    the   recursion  formulas   give  an  exact  rather  than  an 
asymptotic   solution.      This  paper,    therefore,    presents  an  alternate 
approach  to  the   calculation  of  the  low  frequency  fields  propagated  in  the 
terrestrial  waveguide.      Finally  the   solution  obtained  leads  to  an  improved 
method  for  introducing  the  anisotropic  ionosphere   reflection  coefficient 
matrix. 


5- 


2.      Theory  of  Zonal  Harmonic  and  Geometric  Series 

In  the  previous  paper  [johler,    1964]  the  propagated  field  (vertical 
polarization  of  the   source)  was  written  as  the   series, 


Er   =  Er 


,o+I 


'r,    1 


(2.1) 


J=i 


where  the  ground  wave,    Er    0,    can  be  regarded  as  the  zero  order  (j  =  0) 
term  of  the   series, 


[i0c    i0l      f  Tr(„\  r(2)  r(l)n  +  r   ) 


n  =0 


and  a  particular  ionosphere  wave,    j, 


CO 

lipc    i0i    y  r(2)r(l)n  +  r  r 

Er>  J  ~    8n    kp  A    F(n)  Cia  Cla(1  +RJ 


n  =0 


(2.2) 


rr(i)       r(a). 

kliL  .     Zl-6 
L      (3)  (1). 


_Rn] 


j  -i 


T. 


(2.3) 


where  Fn    =  n(n+l)(2n+l)  Pn    (cos    9).      Here  the  order  of  summation  of 
the  zonal  harmonics  and  geometric  series  has  been  interchanged.      This 
does  not  impose  any  serious  mathematical  restrictions,    since  the  j -series 
and  the  n-series  are  absolutely  convergent. 

The  Legendre  function,    Pn(z),    can  be  calculated  exactly  for  integral 
n  from  the   recursion  formula, 


p„  +  1(z)  =  <^i^p„(z)-^ipn.a(z). 


(2.4) 


where, 

Po(z)  =   1  (2.  5) 

and, 

Px(z)   =  z.  (2.6) 

Both  integer  and  non-integer  values  of  n  can  be  used  in  a  formula 
[Bremmer,    1949]  based  on  the   saddle  point  approximation  of  the 
integral  of  Sommerfeld  (letter  z  =  cos   6). 


2  cos 
Pn(cos    0) 


(n+!)9  -  J_j 


V2TT(n+l)   sin   0 


(2.7) 


This  formula  becomes  inaccurate  for  sin   0  ~  0  and  n  small.      This 
formula  (2.  7)  was  checked  with  formulas  (2.  4),    (2.  5)  and  (2.  6)  in  a  com- 
puter subroutine.      The  agreement  was  excellent.      The  recursion  formulas 
were  used,    however,    for  the  computations  illustrated  in  this  paper. 

The  spherical  wave  functions  can  also  be  calculated  from  the 
recursion  formulas, 


(1,2)  2n+j_     (1,3)  r(1»a).    x  ,?    o, 

^n+l       (z)    =    : Cn  (z)    -    Cn-1       (z)>  (2-  8) 


starting  with  the  values, 

Ci1''3  \z)  =   =Fi  exp  [±iz]  (2.9) 

C(_Xi,8)(z)  =    exp  [±iz].  (2.  10) 

Using  a  computation  technique  of  Goldsteing  and  Thayler  [1959],    a 
computer  program  was  written  to  calculate  the  spherical  wave  functions 
(2.  8),    (2.  9),    and  (2.  10).      An  extensive  comparison  was  made  with  a 
program  developed  by  Berry  [  1964a]  using  the  asymptotic  formulas 
(designed  for  the  complex  n-plane)  along  the   real  axis  of  the  n-plane. 


Excellent  agreement  was  found.      The  recursion  formulas  (2.  8),    (2.  9), 
and  (2.  10)  are  of  course  calculated  with  much  greater  speed  and  accuracy 
on  the  computer. 

The  ratio  of  the  derivative  of  a  spherical  wave  function  to  the 
function  can  also  be  determined  by  recursion  [Johler  and  Berry,    1962], 


#'"><.) 

1 

n 

e  ••'(.) " 

n 
z 

(i  ,2  )' 

d.'i  ;  (z) 

Cn-l       (z) 

z 

(2.11) 


where, 


=  ±  i. 


For  complex  argument,    z  =  k2a,    the  Debye  approximation  can  be 
employed, 


(i) 


C;  '  (ksa)  _,     /n(n+l)    _  x 


tfW,        V^2 


(2.12) 


This  very  simple  approximation  becomes  quite  excellent  since  v  n(n+l) 

does  not  grow  as  large  as    | k2 a  | ,    since  the   series  converges  at  n  ~  kx  a. 

is  7(k   a) 
Alternately,    the  function     yn  „  2 — -   can  be   readily  used  in  (2.  11)  by 

Mksa) 

/2\  (2\7 

replacing  the  Cn      and  Q       function  with  \|jn  and  \|^  functions.      Thus, 


CJ1}  (k2a)  K'(k2a) 


(!) 

Ci;(k2a) 


(2.13) 


'n(k2a) 


Equations  (2.  2)  and  (2.  3)  can  be   readily  generalized  for  elevated 
terminals,    with  the   source  S  situated  at  b  =  a  +  hs   above  the  ground  and 
the  observer  O  situated  at  r  =  a  +  h0 ,  instead  of  the  surface  r  =  b  =  a 
illustrated  in  figure   1.      Thus,    the  ground  wave  can  be  written,    r  >  b, 


(1) 


.(i)    ,    A2)    CI  a'   0    1-(a) 


'r,0-     8rT       k2r2b2     A    F(n)   l>lb    +   Gib        (a)   Rv 


n  =0 


CI 


(2.  14) 


and  the  ionosphere  waves  can  be  written, 


E. 


,  1 


8tt 


Iq1       Y    F<ni  [r(l)  +  c!2r(2)  R 

n -0  Gi a 


k^r 


(2) 

x  !cHr.  +^ch)]pjRi_iTi 

Ci  a 


(2.15) 


Here  the  abbreviations,    cfr'3)=  cf  '  2  \^j.  r);   dV  *  *  =  d1 '  *  ^b),  .  .  .  , 
have  been  employed. 

The   series  (2.  2)  is  the  zonal  harmonics  representation  of  the 
ground  wave  and  converges  quite   slowly.      Thus,    approximately  15  kxa 
terms  are  required  to  get  a  reasonably  accurate  answer.      At  10  kc  /  s 
this  would  be  20,  000  terms.      Hence,    the  zonal  harmonic  series  (2.  2) 
should  be  replaced  by  the  classical  ground  wave  theory.      An  efficient 
computer  program  has  been  developed  for  this  purpose  Johler  [1962]. 
Once  the  ground  wave  (2.  2)  has  been  removed  as  a  problem,    the   second 
term  of  (2.  1)  and  the  higher  order  terms  can  be  readily  evaluated  with 
zonal  harmonics.      This  can  be  made  evident  by  introducing  the  Wronskian 

(Johler  and  Berry,    1962  ;     page  768). 


c«l)'     c<2>' 


c(,l)      c<2) 


2i, 


(2.16) 


into  the  factor  1    +  Rn   in  equation  (2.  3).      Thus, 


1  +  R, 


2i 


c(1)c(2) 


■c(2)    k    c(1) 

*sia         f^i_     <s2a 


A2)    '   k2         (!) 


2a 


(2.  17) 


Thus,    the  £^a'  Q{  a    in  (2.  17)  cancels  the  Gla  £la    in  series  (2.  2),    and  as  n 

grows  large,    advantage  cannot  be  taken  of  the  fact  that  Q1&  £la    grows 
large.      This  is  the   reason  for  the  notoriously  slow  convergence  of  the 
harmonic   series   [Johler  and  Berry,     1962],      On  the  other  hand,    the 
ionospheric  waves  converge   rapidly  as  n-4cx  a  since  the  factor  (1  +  Rn) 
replaces  (1  +  Rn)  in  (2.  3).     It  is  noted  that  Rn  -»  -1   as  n  grows  large. 
This  would  lead  to  computation  problems  if  (2.  17)  were  not  employed. 
In  any  case  approximately  1400  terms  at  1 0  kc/s  or  approximately  kx  a 
terms  required.      In  fact,    1500  terms   gives  very  high  accuracy  (>  6 
significant  figures). 

Although  this  is   still  a  large  number  of  terms,    a  simple  recursion 
method  for  calculating  each  term  allows  the  electronic  computer  to  sum 
the   series  as  quickly  as  it  could  calculate  the   smaller  number  of  complex 
terms  used  in  the  evaluation  of  an  integral  in  the   complex  n-plane.      The 
complete  field  can  then  be  calculated  by  a  summation  of  the  geometric 
series  (2.  1)  which  requires  only  a  few  terms  if  the  upper  boundary  is 
lossy.      However,     100  terms  or  more  may  be   required  if  the  boundaries 
are  not  very  lossy,    but  this  adds  little   to  the  computation  time.      In  (2.  3) 
the  factor, 

M)       Jfl) 
Fc(n)  =   -iii-    ^pexp(icp)  (2.18) 

Gla         Gig 
(1)  (2)* 

is  of  unity  magnitude   since  £  *  a'   =  £la      for  real  arguments  kx  a  and  kx  g. 
Also,    as  n  becomes  large,    cp  is  a  phase  angle  which  approaches  zero, 


Lim  Fc(n)  =  1.  (2.  19) 

The  analysis  procedure  described  above  which  in  essence  removes  the 
ground  wave  from  the  problem,    permits  the  practical  use  of  zonal 
harmonics  techniques  to  frequencies  as   great  as  300  kc/s.      Of  course, 
the  linear  dependence  on  kx  a  makes  the  technique  more  efficient  at 
10  kc/s  and  lower.      Below  2  kc/s  the  zonal  harmonics  can  also  be  used 
for  the  ground  wave  without  any  significant  increase  in  computation  time 
on  a  computer.     In  the  case  of  elevated  terminals  (2.  13)  and  (2.  14)  the 

convergence  is  improved  by  the  factor      f  —    J 
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Thus,    the  factor  (see  (2.  14)  and  (2.  15)), 


"sib  +  bi  b  /2\    ^n 

&1  a 


=  [1  +RB]  Cl 


(x) 


(2.  20) 


if  b  =  a  (receiver  only  elevated).      For  n  >  k-j^a,    Johler  and  Berry  [1962] 
used  a  more  quickly  convergent  series  based  on  an  approximation  of 
Watson  [1919],    or  using  (2.  17), 


(2) 


(1  +Rn)(c12r)Rn  +  ^  q[ 


where, 


t1) 


exp  (is 


4ik£a 
2n+l 


exp(icp)  +  R 


^1  r        "sia 


Cl 


(2) 


c(1) 

bl  a 


(2.21) 


Hence,    the  terminals  S  and  O  on  the  ground  require  the  greatest  number 
of  terms.     It  is  also  noted  that  a  large  number  of  cases,    say  100  distances, 
can  be  calculated  simultaneously  since  the  distance  factor  or  Legendre 
function  is  merely  multiplied  by  the  frequency  dependent  factor,    term-by- 
term  to  introduce  the  distance  dependence. 


3.      Plane  Reflection  Coefficients 

The  ionosphere   spherical  reflection  coefficient  (1.  13)  can  be 
replaced  by  the  Fresnel  reflection  coefficient,    which,    for  vertical  polari- 
zation can  be  written, 


cos 


^-(tr^O2 


cos  CPt    + 


1   -(^sincpj 


(3.1) 


where  [Johler  and  Berry,    1964] 


cos  qjj    =  l 


:      Ci 


(")' 


c(*) 


1  n(n+1) 

1     ~   —r. rrr 


(kig)' 


(3.2) 


11- 


or,  

Vn(n+1)  _, 

sin  CD.    ~  -~r* (3.  3) 

1  ki  g 

with  the   stipulation  that  sin  cpt    =  1  if  — — >  1.      This  latter  stipulation 

kig 
causes  only  a  small  loss  in  computation  precision  since  near  n  =  kx  a  the 
series  converges  abruptly  (see  (4.  10)  below),    and  the   spherical  reflection 
coefficient  reaches  a  value  of    -1.      The  plane  reflection  coefficient,    Te  e  , 
exhibits  the   same  asymptotic  behavior  of  n->°°.      The  introduction  of  plane 
Fresnel  reflection  coefficients  permits  use  of  available  reflection 
coefficients  for  an  anisotropic  plane  plasma  developed  by  Johler  and 

Harper  [1962],    Johler  [1963],    and  Wait  and  Walters  [1963].      Fresnel 
reflection  coefficients  are  not  necessary  to  the  isotropic  case  (3.  1) 
except  to  demonstrate  numerically  that  the  procedure  causes  only  small 
loss  in  computation  precision.      The  capability  of  using  Fresnel  reflection 
coefficients  is  an  important  consideration  in  the  anisotropic  ionosphere, 
since  available  reflection  coefficients  are  restricted  to  the  Fresnel  type. 


4.      Inhomogeneous,    Anisotropic  Ionosphere 

The  generalization  of  equation  (2.  3)  to  the  inhomogeneous, 
anisotropic  ionosphere  is   similar  to  methods  introduced  in  previously 
published  work  [Johler  and  Berry,    1964],      A  consequence  of  the  presence 
o f  the  anisotropic  ionosphere  is  the  existence  of  TE  (transverse  electric) 
propagation  in  addition  to  the  TM  (transverse  magnetic)  propagation, 
notwithstanding  the  fact  that  the   source  is  TM.      Equation  (2.  3)  can  be 
re  -written, 


&i  a  &1 


Let  p  =    i  ^ ^LL 

A2)      M) 

tla  tig 


*    =     Re  ,  „     Pn 
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R. 


L  o      -l  J'     n     L 


T  T 

e  e  x  e  m 

T  T 

me  mm 


Pn    =   P 


o    n    7     r  i         oi 
J'J=Lo         iJ 


R 


Then, 


(p^*.^)'"1?^    =  PJ[ 


^CJ      XJ 


y^    zj  J 


(4.2) 


Then  it  can  be  readily  shown  that, 

Cx   =  Te  e 

C2   =  ReT^      +  R   T      T 


o  e       iq       ee       em       me  e       ©  e     ■         m        mm       ein       me 


(4.3) 


These  results  (4.  3)  were  given  by  Johler  [l96l].      These  results  (4.2) 
and  (4.  3)  can  be  deduced  from  equations  given  by  Johler  and  Berry 
[l964]--see  equations  on  page  11;  (43)  and  (52)  on  pages  8  and  9. 

The  solution  for  the  vertical  electric  field,    Er  ,    can  be  written, 
[Berry,    1964b], 


E. 


|i  o  c       l0l 


8TT     k^a4 


CO 

£     F(n)  C^CxVd+Re) 


M+Pn^. 


n  =0 


\1     ~    Me.Jn 


(4.4) 


The  ratio  of  the  determinants  in  (4.  4)  can  be  expanded  in  a  geometric 
type   series  by  Berry  [1964b];  Johler  and  Berry  [l964]--(43)  and  (52), 
pages  8  and  9, 


-13- 


I-Z   +   Pn^l  x 

=    \l-    pA.n^l         U+   Pn^nl  (4-5) 


i-P„^e.n^nl 


or, 


1+    Pn^n 


1    ~    Pn*.      n^n 


=     U    +  (■*+*.,,)     £      (Pn^n^Jpi"1^!.  (4-6) 


5=1 


This,    of  course,    leads  to  (4.  2).      Berry  [1964]  pre -multiplied  (pnTn)j 

by  Bl  ~     which  would  cause  the  second  term  of  C2   to  read  -  Rei?m  Te    Tm  e  . 
Similarly,    the  higher  order  Cj's  would  be  altered.      Whereas  this  permits 
removal  of  a  factor  y .    from  the  integral  (or  in  this  paper  the  summation), 
such  that    Yj   =  Ri_1Cj,    the  procedure,    although  a  good  approximation, 
does  not  seem  to  be  physically  real,    since  the   second  term  is  a  TE  wave 
and  does  not  require  a  TM  reflection  coefficient,    Re  ,    (see  figure  2), 

The  expansion  (4.  6)  converges  absolutely  for  all  n.      Thus  the 
expansion  of  the  denominator  (1.  16)  of  the  left  side  of  (4.  6)  is  the  geo- 
metric series  which  converges  if, 


PlA<nrn|  <  l, 


which  is  certainly  true  for  all  physical  values  of  H e       ,  pn,    and  Tn  ,    if    n 
is  not  too  large.      However, 


Lim    pnTn    =   -   1, 


Lim  7?e  f  n   =   -  1 , 

n  — >  oo 


hence  (4.  5), 


Lim  |i+  PnTn  |    =  0, 


14 


which  multiplies  into  the  geometric  series  (4.  6).      Hence  (4.  6)  converges 
absolutely  for  all  values  of  interest,    i.  e. 


Pn^e 


T 


<   1 


The  geometric  series  representation  permits  the  introduction  of 
local  reflection  coefficients.      This  is  illustrated  by  analogy  to  geometric 
optics  in  figure  1.      Thus,    the  first  ionosphere  wave,    j  =  1,    has  a 
reflection  point  [l,  1  ],    located  in  the  middle  of  the  propagation  path.     In 
the  region  of  this  point,    it  can  be  conjectured,    most  of  the  reflection 
occurs.      This  becomes  evident  from  an  examination  of  the  angle  of 
incidence,    cpt  ,    on  the  ionosphere  of  the   spherical  wave  [Johler  and  Berry, 
1964],    determined  by 


cos  cpj    r  Rei 


or 


r(2)' 


Re 


n(n+  1) 
(klg)2 


(4.7) 


sm 


Vn(n+1) 


For  cpt    ~  0,    or  n  small,    the  contribution  of  the  terms  of  the  harmonic 
series  is  small.      As  sin  cpt    approaches  unity,    the  field  is  determined. 
Thereafter,    the   series  converges  rapidly  as  a  consequence  of  the  behavior 
of  (2.  17)  previously  discussed.      Therefore,    within  a  small  range  of 
values  near  kx  g  the  field  is  primarily  determined.      Similarly,    the  local 
reflection  points  [2,  1],    [2,2],    [3,1],    [3,2],    [3,3],    [4,1],    [4,2], 
[4,  3],    and  [4,  4]  can  be  defined.      The  anisotropy  and  non-homogeneity, 
of  course,    causes  the  calculation  of  Cj    to  be  quite  complicated.      Thus, 
instead  of  (4.  2),    use, 


11       (Pn.K^.K^.n.KKPn,^  )=      fl        PK     [^        *< 

=2  l.  y  i         ^  i 


K  =2 


K  =1 


(4.8) 


where  for  j   =  1 , 


(Pn.l^n,!)    =   Pi 


Cl       Xll 
y,        ZXJ 
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Explicit  expression  has  been  given  by  Johler  [l96l]  which  can  be  written 

Ca   =  Tee(l,l) 

C2   =  Tee(2,l)Tee(2,2)Re(2,l)  +  Rn  (2,  1 )  Te  ,  (2,  1 )  Tm  e  (2,  2) 

C3  =  Re  (3,  2)  Rffl  (3,  1)  Tee  (3,  3)  Tem  (3,  1)  T,  .  (3,  2) 

+  Re(3,l)Rm(3,2)Tee(3)l)TeB(3,2)TIBe(3,3) 

+  Re(3,l)Re(3,2)Tee(3,l)Tee(3,2)Tee(3,3) 

+  Rffl  (3,  1)  RB  (3,  2)  Tmffl(3,  2)  Tem  (3,  1) Tffle(3,  3)  (4.  9) 

Since  higher  order  C.'s  are  quite  complicated,    it  is  perhaps   simpler 
to  calculate  C,   from  the  general  matrix  formula  (4.  8).      The  detailed 
reflection,    figure   1,    can  also  be  altered  by  the  introduction  of  a  local 
reflection  height  at  each  point  [1,1],    [2,  1  ],    and  [2,  2]  •  •  •  .      This  permits 
an  extension  of  the  theory  to  take  account  of  irregularities  of  the  iono- 
sphere.     The  limitations  of  this  procedure  are  not  treated  in  this  paper, 
but  the  procedure  is  intuitively  evident  from  the   rays  depicted  in  figure  1. 
It  should  be  emphasized  here  that  the  analysis  presented  here  does  not 
depend  upon  the  geometric -optics  depicted  for  tutorial  purposes  in  figure   1. 
The  solution  is  exact  since  the  waves  fill  all  of  the   space   shown.      The 
significance  of  the  geometric -optical  ray  limit  has  been  discussed  else- 
where,   for  example,    Johler  [1964]  and  Wait  [1961]. 

The  error  introduced  in  (4.  9)  by  using  available  plane  (Fresnel) 
reflection  coefficients  for  real  angles  of  incidence,    cpt  ,    is  small  compared 
with  errors  introduced  by  lack  of  knowledge  of  the  ionosphere  model,    for 
example,    but  the  error  is  not  mathematically  negligible.      Some  improve- 
ment can  be  obtained  by  multiplication  of  (4.  8)  by  a  spherical  correction 
factor,    CF,    based  on  the  spherical  reflection  coefficient  for  infinitely 
conducting  ionosphere, 


These  functions  are,    of  course,    available  from  the  computer  code  written 
for  (2.  2),    and  (2.  3).      However,    the  best  procedure  would  be  the  introduc- 
tion of  the  complex  angle  of  incidence,    cpt  f    (see  (3.  2))    into  the  calculation 
of  the   Fresnel  ionosphere  reflection  coefficients  directly.      Thus,    cos  cpt 
in  the  analysis  procedure  of  Johler  and  Harper  [1962]  is  replaced  by  the 
complex  number  after  equality  sign  in  (3.  2). 
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5.      Discussion  of  a  Sample  Calculation 

Figures  2  and  3  illustrate  a  sample  calculation  using  the  zonal 

harmonic- -geometric  series  double  summation  at  1 0  kc/s  and  20  kc/s. 

An  entire  calculation  of  this  type  requires  approximately  two  minutes  at 

10  kc/s,    (or  perhaps  less  with  optimization  of  the  computer  code)  on  the 

electronic  computer.     The  ground  wave,    j   =  0,    was  calculated  with  the 

aid  of  a  standard  computer  code  using  the  classical  ground  wave  theory. 

At  frequencies  less  Lhan  2  kc/s,    (2.  2)  becomes  efficient,    and  should  be 

employed.      The  final  answer,    E,    employed  20  terms.     However,    graphi- 

i 
cal  accuracy  can  be  obtained  with  fewer  terms.      Since  kx  g  ~  1340,    1400 

terms  of  the  harmonic  series  were  employed  in  figure  2  and  2800  in 

1400 

figure  3.     In  practice,    it  is  usually  only  necessary  to  form      E        at  1  0 

28(30  n=10QQ 

kc/s  and        E         at  20  kc/s.  .  .  ,    (2.  3)  since  the   small  n  terms  ordinarily 

n=20Q0 

contribute  a  negligible  amount  to  the  total  field.      This  results  in  approxi- 
mately 1/5  reduction  in  computer  time.      It  is,    however,    necessary  to 
recur se  from  n  =  0  in  (2.  8)  and  (2.  11)  but  (2.  7)  can  in  this  case  be  used 
to  replace  (2.  4).      The  E  curve  and  curves  j   =  0,    1,    2,    3  agree  with 

previous  results  of  Berry  [  1964b]  using  the  integral  form.      The  number 
of  distances  calculated  by  any  summation  can  be  increased  to  100  or  more 
without  appreciably  increasing  the  computation  time,    since  Pn(cos   0)  is 
a  distance  multiplier  for  each  frequency  dependent  term  which  can  be 
calculated  with  great  speed  using  (2.  4)  or  (2.  7).      The  latter  is  quite 
accurate,    except  at  extremely  low  frequencies  where   small  n  is  impor- 
tant or  at  distances,    sin  6  ~  0. 

It  is  of  interest  to  note  the   characteristic  behavior  of  each  term  of 
the  geometric  series  at  great  distance.      Thus,    the  attenuation  with  dis- 
tance approaches  a  ground  wave  slope.      However,    since  full  wave  type 
solution  is  employed,    this  attenuation  with  distance  does  differ  slightly 
from  the  ground  wave  curve. 

The  complicated  behavior  of  the  particular  case,    f  =  10  kc/s, 
N  -  56  el/cc,    j  =  4,    5,    6,    in  the  region  immediately  preceding  the 
diffraction  region  is  of  particular  interest.      The  apparent  standing  wave 
appears  to  be  a  consequence  of  the  fact  that  a  full  wave  solution  has  been 
employed.      Thus,    the  geometric  optics  solution  would  yield  a  simple 
pseudo-Brewster  angle  type  cusp  in  the   curve.      Here,    however,    a  multi- 
plicity of  waves  exist  and  would  be  expected  under   some   circumstance  to 
yield  a  more  complicated  behavior  near  the  pseudo-Brewster  angle. 
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6.      Conclusions 

A  new  approach  to  the  calculation  of  LF,  VLF,  ELF  fields  has 
been  demonstrated.      The   solution  permits  introduction  of  existing  plane, 
anisotropic  reflection  coefficients  locally  as  a  function  of  a  real  or  com- 
plex angle  of  incidence.      This  permits  retention  of  the  more  exact  method 
of  calculation  in  which  the  anisotropic  ionosphere  reflection  coefficient 
matrix  is  left  inside  the  summation  or  integration  process  without  an 
excessive  increase  in  computer  time.      The  method  appears  to  be  tractable 
up  to  frequencies  of  approximately  300  kc/s  although  the  efficiency  is 
reduced  at  LF  and  MF.      The  method  is  very  efficient  at  VLF  and  ELF. 
The  method  also  can  be  developed  to  introduce  local  irregularities  and 
local  electrical  constant  changes  at  both  the  ground  and  ionosphere. 


The  author  is  indebted  to  his  colleague,    R.    L.    Lewis,    for  the 
design  of  the  computer  subroutines  based  on  the  recursion  formulas. 
The  author  is  also  indebted  to  L.    A.    Berry  for  the  design  of  the 
asymptotic  subroutines  for  the  Hankel  functions  against  which  the 
recursion  formulas  were  checked.      Finally,    the  author  wishes  to  note 
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zonal  harmonic  program  design. 
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Figure  2. 


Illustrating  components  of  the  rigorous  geometric  series 

j   =  0,     1,    2,    3.  .  .    together  with  complete  field,    2,    for  terres 

trial  waveguide  propagation  at  10  kc/s.  i 
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Figure  3.       Illustrating  components  of  the  rigorous  geometric  series 


j  =  0,    1,    2,    3. 


together  with  complete  field,    Z,    for  terres 


trial  waveguide  propagation  at  20  kc/s. 
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